Write a custom CUDA kernel to optimize `Generalized Dice Loss` (GDL).

Formula:
GDL = 1 - 2 * (Sum_l w_l * Sum_n (r_ln * p_ln)) / (Sum_l w_l * Sum_n (r_ln + p_ln))
Where:
- l is class index, n is spatial index (pixels/voxels).
- p_ln is softmax probability.
- r_ln is one-hot ground truth.
- w_l = 1 / (Sum_n r_ln)^2.

Problem Analysis:
1. Memory Bandwidth: Standard implementation involves Softmax, One-hot generation, Summation per class for weights, and then weighted intersection/union sums. This requires multiple passes over the large (N, C, Spatial) tensors.
2. Intermediate Storage: Storing probability maps and one-hot targets consumes significant memory.

Optimization Strategy: Fused Softmax-Reduction Kernel

1. One-Block-per-Sample: Each block processes one image/volume in the batch.

2. On-the-fly Calculation:
   - Compute Softmax probabilities `p` from logits on-the-fly using cached Max/SumExp.
   - Generate `r` (one-hot) from target indices on-the-fly.

3. Fused Accumulation:
   - Iterate over all spatial pixels `n` and classes `l`.
   - Maintain 3 accumulators per class in Shared Memory/Registers:
     - `sum_r`: sum(r_ln) [for weights]
     - `sum_inter`: sum(r_ln * p_ln)
     - `sum_union`: sum(r_ln + p_ln)
   
   Wait, `w_l` depends on `sum_r` across the WHOLE spatial dimension. So we must accumulate `sum_r` for all pixels first? 
   Yes, standard GDL weights depend on the GT volume.
   However, `sum_inter` and `sum_union` also require summation over `n`.
   We can accumulate `sum_r`, `sum_inter`, `sum_union` simultaneously in one pass over spatial dimensions.

4. Shared Memory Reduction:
   - Use atomicAdd or tree reduction in Shared Memory to aggregate these sums for each class across threads.

5. Final Composition:
   - Thread 0 reads the aggregated per-class sums.
   - Computes `w_l = 1 / (sum_r^2 + eps)`.
   - Computes numerator `2 * sum(w_l * sum_inter)` and denominator `sum(w_l * sum_union)`.
   - Writes the final loss per sample.
  
Here's an example to show you the syntax of inline embedding custom CUDA operators in torch: The example given architecture is:   
  
```python
import torch
import torch.nn as nn
import torch.nn.functional as F

BATCH_SIZE = 16
NUM_CLASSES = 4
DEPTH = 32
HEIGHT = 128
WIDTH = 128
SPATIAL_DIM = DEPTH * HEIGHT * WIDTH 
SHAPE_LOGITS = (BATCH_SIZE, NUM_CLASSES, DEPTH, HEIGHT, WIDTH)
SHAPE_TARGET = (BATCH_SIZE, DEPTH, HEIGHT, WIDTH)

EPS = 1e-6
REDUCTION = 'none'

class GeneralizedDiceLoss(nn.Module):
    """
    Generalized Dice Loss (Sudre et al. MICCAI 2017)
    https://arxiv.org/pdf/1707.03237
    """
    def __init__(self, eps=1e-6, reduction='mean'):
        super(GeneralizedDiceLoss, self).__init__()
        self.eps = eps
        self.reduction = reduction

    def forward(self, logits: torch.Tensor, targets: torch.Tensor) -> torch.Tensor:
        # logits: (B, C, D, H, W)
        # targets: (B, D, H, W) -> indices
        probs = F.softmax(logits, dim=1)
        
        # One-hot encoding
        targets_onehot = F.one_hot(targets.long(), num_classes=logits.shape[1])
        targets_onehot = targets_onehot.permute(0, 4, 1, 2, 3).float()
        
        # Compute weights w_l = 1 / (sum(r_ln)^2)
        sum_r = targets_onehot.sum(dim=(2, 3, 4)) # (B, C)
        weights = 1.0 / (sum_r * sum_r + self.eps)
        
        # Compute Intersection & Union
        # Intersection: r * p
        intersection = (targets_onehot * probs).sum(dim=(2, 3, 4)) # (B, C)
        # Union: r + p
        union = (targets_onehot + probs).sum(dim=(2, 3, 4)) # (B, C)
        
        # Weighted Sum
        # Numerator: 2 * sum_l (w_l * inter_l)
        # Denominator: sum_l (w_l * union_l)
        
        numerator = 2.0 * (weights * intersection).sum(dim=1)
        denominator = (weights * union).sum(dim=1)
        
        dice_score = numerator / (denominator + self.eps)
        loss = 1.0 - dice_score
        
        if self.reduction == 'mean':
            return loss.mean()
        elif self.reduction == 'sum':
            return loss.sum()
        return loss

class Model(nn.Module):
    def __init__(self, eps=1e-6, reduction='none'):
        super(Model, self).__init__()
        self.loss_fn = GeneralizedDiceLoss(eps=eps, reduction=reduction)
    
    def forward(self, logits, targets):
        return self.loss_fn(logits, targets)

def get_inputs():
    logits = torch.randn(SHAPE_LOGITS, dtype=torch.float32)
    targets = torch.randint(0, NUM_CLASSES, SHAPE_TARGET, dtype=torch.long)
    return [logits.contiguous(), targets.contiguous()]

def get_init_inputs():
    return [EPS, REDUCTION]